{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# scikit-image advanced panorama tutorial\n",
    "\n",
    "Enhanced from the original demo as featured in [the scikit-image paper](https://peerj.com/articles/453/).\n",
    "\n",
    "Multiple overlapping images of the same scene, combined into a single image, can yield amazing results. This tutorial will illustrate how to accomplish panorama stitching using scikit-image, from loading the images to cleverly stitching them together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First things first\n",
    "\n",
    "Import NumPy and matplotlib, then define a utility function to compare multiple images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def compare(*images, **kwargs):\n",
    "    \"\"\"\n",
    "    Utility function to display images side by side.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    image0, image1, image2, ... : ndarrray\n",
    "        Images to display.\n",
    "    labels : list\n",
    "        Labels for the different images.\n",
    "    \"\"\"\n",
    "    f, axes = plt.subplots(1, len(images), **kwargs)\n",
    "    axes = np.array(axes, ndmin=1)\n",
    "    \n",
    "    labels = kwargs.pop('labels', None)\n",
    "    if labels is None:\n",
    "        labels = [''] * len(images)\n",
    "    \n",
    "    for n, (image, label) in enumerate(zip(images, labels)):\n",
    "        axes[n].imshow(image, interpolation='nearest', cmap='gray')\n",
    "        axes[n].set_title(label)\n",
    "        axes[n].axis('off')\n",
    "    \n",
    "    f.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data\n",
    "\n",
    "The ``ImageCollection`` class provides an easy and efficient way to load and represent multiple images. Images in the ``ImageCollection`` are not only read from disk when accessed.\n",
    "\n",
    "Load a series of images into an ``ImageCollection`` with a wildcard, as they share similar names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage.io as io\n",
    "\n",
    "pano_imgs = io.ImageCollection('../../images/pano/JDW_03*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect these images using the convenience function `compare()` defined earlier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare(*pano_imgs, figsize=(12, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Credit: Images of Private Arch and the trail to Delicate Arch in Arches National Park, USA, taken by Joshua D. Warner.<br>\n",
    "License: CC-BY 4.0\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0. Pre-processing\n",
    "\n",
    "This stage usually involves one or more of the following:\n",
    "* Resizing, often downscaling with fixed aspect ratio\n",
    "* Conversion to grayscale, as many feature descriptors are not defined for color images\n",
    "* Cropping to region(s) of interest\n",
    "\n",
    "For convenience our example data is already resized smaller, and we won't bother cropping. However, they are presently in color so coversion to grayscale with `skimage.color.rgb2gray` is appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.color import rgb2gray\n",
    "\n",
    "pano0, pano1, pano2 = [rgb2gray(im) for im in pano_imgs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the results\n",
    "compare(pano0, pano1, pano2, figsize=(12, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Feature detection and matching\n",
    "\n",
    "We need to estimate a projective transformation that relates these images together. The steps will be\n",
    "\n",
    "1. Define one image as a _target_ or _destination_ image, which will remain anchored while the others are warped\n",
    "2. Detect features in all three images\n",
    "3. Match features from left and right images against the features in the center, anchored image.\n",
    "\n",
    "In this three-shot series, the middle image `pano1` is the logical anchor point.\n",
    "\n",
    "We detect \"Oriented FAST and rotated BRIEF\" (ORB) features in both images. \n",
    "\n",
    "**Note:** For efficiency, in this tutorial we're finding 800 keypoints. The results are good but small variations are expected. If you need a more robust estimate in practice, run multiple times and pick the best result _or_ generate additional keypoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import ORB\n",
    "\n",
    "# Initialize ORB\n",
    "# 800 keypoints is large enough for robust results, \n",
    "# but low enough to run within a few seconds. \n",
    "orb = ORB(n_keypoints=800, fast_threshold=0.05)\n",
    "\n",
    "# Detect keypoints in pano0\n",
    "orb.detect_and_extract(pano0)\n",
    "keypoints0 = orb.keypoints\n",
    "descriptors0 = orb.descriptors\n",
    "\n",
    "# Detect keypoints in pano1\n",
    "orb.detect_and_extract(pano1)\n",
    "keypoints1 = orb.keypoints\n",
    "descriptors1 = orb.descriptors\n",
    "\n",
    "# Detect keypoints in pano2\n",
    "orb.detect_and_extract(pano2)\n",
    "keypoints2 = orb.keypoints\n",
    "descriptors2 = orb.descriptors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Match features from images 0 <-> 1 and 1 <-> 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import match_descriptors\n",
    "\n",
    "# Match descriptors between left/right images and the center\n",
    "matches01 = match_descriptors(descriptors0, descriptors1, cross_check=True)\n",
    "matches12 = match_descriptors(descriptors1, descriptors2, cross_check=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect these matched features side-by-side using the convenience function ``skimage.feature.plot_matches``. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import plot_matches\n",
    "fig, ax = plt.subplots(1, 1, figsize=(15, 12))\n",
    "\n",
    "# Best match subset for pano0 -> pano1\n",
    "plot_matches(ax, pano0, pano1, keypoints0, keypoints1, matches01)\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of these line up similarly, but it isn't perfect. There are a number of obvious outliers or false matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 12))\n",
    "\n",
    "# Best match subset for pano2 -> pano1\n",
    "plot_matches(ax, pano1, pano2, keypoints1, keypoints2, matches12)\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to above, decent signal but numerous false matches.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Transform estimation\n",
    "\n",
    "To filter out the false matches, we apply RANdom SAmple Consensus (RANSAC), a powerful method of rejecting outliers available in ``skimage.transform.ransac``.  The transformation is estimated iteratively, based on randomly chosen subsets, finally selecting the model which corresponds best with the majority of matches.\n",
    "\n",
    "We need to do this twice, once each for the transforms left -> center and right -> center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import ProjectiveTransform\n",
    "from skimage.measure import ransac\n",
    "\n",
    "# Select keypoints from \n",
    "#   * source (image to be registered): pano0\n",
    "#   * target (reference image): pano1, our middle frame registration target\n",
    "src = keypoints0[matches01[:, 0]][:, ::-1]\n",
    "dst = keypoints1[matches01[:, 1]][:, ::-1]\n",
    "\n",
    "model_robust01, inliers01 = ransac((src, dst), ProjectiveTransform,\n",
    "                                   min_samples=4, residual_threshold=1, max_trials=300)\n",
    "\n",
    "# Select keypoints from \n",
    "#   * source (image to be registered): pano2\n",
    "#   * target (reference image): pano1, our middle frame registration target\n",
    "src = keypoints2[matches12[:, 1]][:, ::-1]\n",
    "dst = keypoints1[matches12[:, 0]][:, ::-1]\n",
    "\n",
    "model_robust12, inliers12 = ransac((src, dst), ProjectiveTransform,\n",
    "                                   min_samples=4, residual_threshold=1, max_trials=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `inliers` returned from RANSAC select the best subset of matches. How do they look?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 12))\n",
    "\n",
    "# Best match subset for pano0 -> pano1\n",
    "plot_matches(ax, pano0, pano1, keypoints0, keypoints1, matches01[inliers01])\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(15, 12))\n",
    "\n",
    "# Best match subset for pano2 -> pano1\n",
    "plot_matches(ax, pano1, pano2, keypoints1, keypoints2, matches12[inliers12])\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the false matches are rejected!\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Warping\n",
    "\n",
    "Next, we produce the panorama itself. We must _warp_, or transform, two of the three images so they will properly align with the stationary image.\n",
    "\n",
    "### Extent of output image\n",
    "The first step is to find the shape of the output image to contain all three transformed images. To do this we consider the extents of all warped images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import SimilarityTransform\n",
    "\n",
    "# Shape of middle image, our registration target\n",
    "r, c = pano1.shape[:2]\n",
    "\n",
    "# Note that transformations take coordinates in (x, y) format,\n",
    "# not (row, column), in order to be consistent with most literature\n",
    "corners = np.array([[0, 0],\n",
    "                    [0, r],\n",
    "                    [c, 0],\n",
    "                    [c, r]])\n",
    "\n",
    "# Warp the image corners to their new positions\n",
    "warped_corners01 = model_robust01(corners)\n",
    "warped_corners12 = model_robust12(corners)\n",
    "\n",
    "# Find the extents of both the reference image and the warped\n",
    "# target image\n",
    "all_corners = np.vstack((warped_corners01, warped_corners12, corners))\n",
    "\n",
    "# The overall output shape will be max - min\n",
    "corner_min = np.min(all_corners, axis=0)\n",
    "corner_max = np.max(all_corners, axis=0)\n",
    "output_shape = (corner_max - corner_min)\n",
    "\n",
    "# Ensure integer shape with np.ceil and dtype conversion\n",
    "output_shape = np.ceil(output_shape[::-1]).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply estimated transforms\n",
    "\n",
    "Warp the images with `skimage.transform.warp` according to the estimated models. A shift, or _translation_ is needed to place as our middle image in the middle - it isn't truly stationary.\n",
    "\n",
    "Values outside the input images are initially set to -1 to distinguish the \"background\", which is identified for later use.\n",
    "\n",
    "**Note:** ``warp`` takes the _inverse_ mapping as an input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import warp\n",
    "\n",
    "# This in-plane offset is the only necessary transformation for the middle image\n",
    "offset1 = SimilarityTransform(translation= -corner_min)\n",
    "\n",
    "# Translate pano1 into place\n",
    "pano1_warped = warp(pano1, offset1.inverse, order=3,\n",
    "                    output_shape=output_shape, cval=-1)\n",
    "\n",
    "# Acquire the image mask for later use\n",
    "pano1_mask = (pano1_warped != -1)  # Mask == 1 inside image\n",
    "pano1_warped[~pano1_mask] = 0      # Return background values to 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warp left panel into place"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Warp pano0 (left) to pano1\n",
    "transform01 = (model_robust01 + offset1).inverse\n",
    "pano0_warped = warp(pano0, transform01, order=3,\n",
    "                    output_shape=output_shape, cval=-1)\n",
    "\n",
    "pano0_mask = (pano0_warped != -1)  # Mask == 1 inside image\n",
    "pano0_warped[~pano0_mask] = 0      # Return background values to 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warp right panel into place"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Warp pano2 (right) to pano1 \n",
    "transform12 = (model_robust12 + offset1).inverse\n",
    "pano2_warped = warp(pano2, transform12, order=3,\n",
    "                    output_shape=output_shape, cval=-1)\n",
    "\n",
    "pano2_mask = (pano2_warped != -1)  # Mask == 1 inside image\n",
    "pano2_warped[~pano2_mask] = 0      # Return background values to 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the warped images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare(pano0_warped, pano1_warped, pano2_warped, figsize=(12, 10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Combining images the easy (and bad) way\n",
    "\n",
    "This method simply \n",
    "\n",
    "1. sums the warped images\n",
    "2. tracks how many images overlapped to create each  point\n",
    "3. normalizes the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the three images together. This could create dtype overflows!\n",
    "# We know they are are floating point images after warping, so it's OK.\n",
    "merged = (pano0_warped + pano1_warped + pano2_warped)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Track the overlap by adding the masks together\n",
    "overlap = (pano0_mask * 1.0 +  # Multiply by 1.0 for bool -> float conversion\n",
    "           pano1_mask + \n",
    "           pano2_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize through division by `overlap` - but ensure the minimum is 1\n",
    "normalized = merged / np.maximum(overlap, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, view the results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "ax.imshow(normalized, cmap='gray')\n",
    "\n",
    "plt.tight_layout()\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "<div style=\"height: 400px;\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What happened?!** Why are there nasty dark lines at boundaries, and why does the middle look so blurry?\n",
    "\n",
    "\n",
    "The **lines are artifacts (boundary effect) from the warping method**. When the image is warped with interpolation, edge pixels containing part image and part background combine these values. We would have bright lines if we'd chosen `cval=2` in the `warp` calls (try it!), but regardless of choice there will always be discontinuities.\n",
    "\n",
    "...Unless you use `order=0` in `warp`, which is nearest neighbor. Then edges are perfect (try it!). But who wants to be limited to an inferior interpolation method? \n",
    "\n",
    "Even then, it's blurry! Is there a better way?\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Stitching images along a minimum-cost path\n",
    "\n",
    "Let's step back a moment and consider: Is it even reasonable to blend pixels?\n",
    "\n",
    "Take a look at a difference image, which is just one image subtracted from the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,12))\n",
    "\n",
    "# Generate difference image and inspect it\n",
    "difference_image = pano0_warped - pano1_warped\n",
    "ax.imshow(difference_image, cmap='gray')\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The surrounding flat gray is zero. _A perfect overlap would show no structure!_ \n",
    "\n",
    "Instead, the overlap region matches fairly well in the middle... but off to the sides where things start to look a little embossed, a simple average blurs the result. This caused the blurring in the previous, method (look again). _Unfortunately, this is almost always the case for panoramas!_\n",
    "\n",
    "How can we fix this?\n",
    "\n",
    "Let's attempt to find a vertical path through this difference image which stays as close to zero as possible. If we use that to build a mask, defining a transition between images, the result should appear _seamless_.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Seamless image stitching with Minimum-Cost Paths and `skimage.graph`\n",
    "\n",
    "Among other things, `skimage.graph` allows you to\n",
    "* start at any point on an array\n",
    "* find the path to any other point in the array\n",
    "* the path found _minimizes_ the sum of values on the path.\n",
    "\n",
    "\n",
    "The array is called a _cost array_, while the path found is a _minimum-cost path_ or **MCP**.\n",
    "\n",
    "To accomplish this we need\n",
    "\n",
    "* Starting and ending points for the path\n",
    "* A cost array (a modified difference image)\n",
    "\n",
    "This method is so powerful that, with a carefully constructed cost array, the seed points are essentially irrelevant. It just works!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define seed points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmax = output_shape[0] - 1\n",
    "cmax = output_shape[1] - 1\n",
    "\n",
    "# Start anywhere along the top and bottom, left of center.\n",
    "mask_pts01 = [[0,    cmax // 3],\n",
    "              [rmax, cmax // 3]]\n",
    "\n",
    "# Start anywhere along the top and bottom, right of center.\n",
    "mask_pts12 = [[0,    2*cmax // 3],\n",
    "              [rmax, 2*cmax // 3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct cost array\n",
    "\n",
    "This utility function exists to give a \"cost break\" for paths from the edge to the overlap region.\n",
    "\n",
    "We will visually explore the results shortly. Examine the code later - for now, just use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import flood_fill\n",
    "\n",
    "def generate_costs(diff_image, mask, vertical=True, gradient_cutoff=2.):\n",
    "    \"\"\"\n",
    "    Ensures equal-cost paths from edges to region of interest.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    diff_image : ndarray of floats\n",
    "        Difference of two overlapping images.\n",
    "    mask : ndarray of bools\n",
    "        Mask representing the region of interest in ``diff_image``.\n",
    "    vertical : bool\n",
    "        Control operation orientation.\n",
    "    gradient_cutoff : float\n",
    "        Controls how far out of parallel lines can be to edges before\n",
    "        correction is terminated. The default (2.) is good for most cases.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    costs_arr : ndarray of floats\n",
    "        Adjusted costs array, ready for use.\n",
    "    \"\"\"\n",
    "    if vertical is not True:\n",
    "        return tweak_costs(diff_image.T, mask.T, vertical=vertical,\n",
    "                           gradient_cutoff=gradient_cutoff).T\n",
    "\n",
    "    # Start with a high-cost array of 1's\n",
    "    costs_arr = np.ones_like(diff_image)\n",
    "\n",
    "    # Obtain extent of overlap\n",
    "    row, col = mask.nonzero()\n",
    "    cmin = col.min()\n",
    "    cmax = col.max()\n",
    "    shape = mask.shape\n",
    "\n",
    "    # Label discrete regions\n",
    "    labels = mask.copy().astype(np.uint8)\n",
    "    cslice = slice(cmin, cmax + 1)\n",
    "    submask = np.ascontiguousarray(labels[:, cslice])\n",
    "    submask = flood_fill(submask, (0, 0), 2)\n",
    "    submask = flood_fill(submask, (shape[0]-1, 0), 3)\n",
    "    labels[:, cslice] = submask\n",
    "\n",
    "    # Find distance from edge to region\n",
    "    upper = (labels == 2).sum(axis=0).astype(np.float64)\n",
    "    lower = (labels == 3).sum(axis=0).astype(np.float64)\n",
    "\n",
    "    # Reject areas of high change\n",
    "    ugood = np.abs(np.gradient(upper[cslice])) < gradient_cutoff\n",
    "    lgood = np.abs(np.gradient(lower[cslice])) < gradient_cutoff\n",
    "\n",
    "    # Give areas slightly farther from edge a cost break\n",
    "    costs_upper = np.ones_like(upper)\n",
    "    costs_lower = np.ones_like(lower)\n",
    "    costs_upper[cslice][ugood] = upper[cslice].min() / np.maximum(upper[cslice][ugood], 1)\n",
    "    costs_lower[cslice][lgood] = lower[cslice].min() / np.maximum(lower[cslice][lgood], 1)\n",
    "\n",
    "    # Expand from 1d back to 2d\n",
    "    vdist = mask.shape[0]\n",
    "    costs_upper = costs_upper[np.newaxis, :].repeat(vdist, axis=0)\n",
    "    costs_lower = costs_lower[np.newaxis, :].repeat(vdist, axis=0)\n",
    "\n",
    "    # Place these in output array\n",
    "    costs_arr[:, cslice] = costs_upper[:, cslice] * (labels[:, cslice] == 2)\n",
    "    costs_arr[:, cslice] +=  costs_lower[:, cslice] * (labels[:, cslice] == 3)\n",
    "\n",
    "    # Finally, place the difference image\n",
    "    costs_arr[mask] = diff_image[mask]\n",
    "\n",
    "    return costs_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use this function to generate the cost array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with the absolute value of the difference image.\n",
    "# np.abs is necessary because we don't want negative costs!\n",
    "costs01 = generate_costs(np.abs(pano0_warped - pano1_warped),\n",
    "                         pano0_mask & pano1_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allow the path to \"slide\" along top and bottom edges to the optimal horizontal position by setting top and bottom edges to zero cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "costs01[0,  :] = 0\n",
    "costs01[-1, :] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our cost array now looks like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 12))\n",
    "\n",
    "ax.imshow(costs01, cmap='gray', interpolation='none')\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tweak we made with `generate_costs` is subtle but important. Can you see it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the minimum-cost path (MCP)\n",
    "\n",
    "Use `skimage.graph.route_through_array` to find an optimal path through the cost array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.graph import route_through_array\n",
    "\n",
    "# Arguments are:\n",
    "#   cost array\n",
    "#   start pt\n",
    "#   end pt\n",
    "#   can it traverse diagonally\n",
    "pts, _ = route_through_array(costs01, mask_pts01[0], mask_pts01[1], fully_connected=True)\n",
    "\n",
    "# Convert list of lists to 2d coordinate array for easier indexing\n",
    "pts = np.array(pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did it work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "# Plot the difference image\n",
    "ax.imshow(pano0_warped - pano1_warped, cmap='gray')\n",
    "\n",
    "# Overlay the minimum-cost path\n",
    "ax.plot(pts[:, 1], pts[:, 0])  \n",
    "\n",
    "plt.tight_layout()\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks like a great seam to stitch these images together - the path looks very close to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Irregularities\n",
    "\n",
    "Due to the random element in the RANSAC transform estimation, everyone will have a slightly different path. **Your path will look different from mine, and different from your neighbor's.** That's expected! _The awesome thing about MCP is that everyone just calculated the best possible path to stitch together their unique transforms!_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filling the mask\n",
    "\n",
    "Turn that path into a mask, which will be 1 where we want the left image to show through and zero elsewhere. We need to fill the left side of the mask with ones over to our path.\n",
    "\n",
    "**Note**: This is the inverse of NumPy masked array conventions (``numpy.ma``), which specify a negative mask (mask == bad/missing) rather than a positive mask as used here (mask == good/selected).\n",
    "\n",
    "Place the path into a new, empty array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with an array of zeros and place the path\n",
    "mask0 = np.zeros_like(pano0_warped, dtype=np.uint8)\n",
    "mask0[pts[:, 0], pts[:, 1]] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ensure the path appears as expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(11, 11))\n",
    "\n",
    "# View the path in black and white\n",
    "ax.imshow(mask0, cmap='gray')\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Label the various contiguous regions in the image using `skimage.measure.label`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import flood_fill\n",
    "\n",
    "# Labeling starts with one at point (0, 0)\n",
    "mask0 = flood_fill(mask0, (0, 0), 1, connectivity=1)\n",
    "\n",
    "# The result\n",
    "plt.imshow(mask0, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks great!\n",
    "\n",
    "### Rinse and repeat\n",
    "\n",
    "Apply the same principles to images 1 and 2: first, build the cost array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with the absolute value of the difference image.\n",
    "# np.abs necessary because we don't want negative costs!\n",
    "costs12 = generate_costs(np.abs(pano1_warped - pano2_warped),\n",
    "                         pano1_mask & pano2_mask)\n",
    "\n",
    "# Allow the path to \"slide\" along top and bottom edges to the optimal \n",
    "# horizontal position by setting top and bottom edges to zero cost\n",
    "costs12[0,  :] = 0\n",
    "costs12[-1, :] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Add an additional constraint this time**, to prevent this path crossing the prior one!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "costs12[mask0 > 0] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "ax.imshow(costs12, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your results may look slightly different.\n",
    "\n",
    "Compute the minimal cost path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arguments are:\n",
    "#   cost array\n",
    "#   start pt\n",
    "#   end pt\n",
    "#   can it traverse diagonally\n",
    "pts, _ = route_through_array(costs12, mask_pts12[0], mask_pts12[1], fully_connected=True)\n",
    "\n",
    "# Convert list of lists to 2d coordinate array for easier indexing\n",
    "pts = np.array(pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify a reasonable result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "# Plot the difference image\n",
    "ax.imshow(pano1_warped - pano2_warped, cmap='gray')\n",
    "\n",
    "# Overlay the minimum-cost path\n",
    "ax.plot(pts[:, 1], pts[:, 0]);\n",
    "\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the mask by placing the path in a new array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask2 = np.zeros_like(pano0_warped, dtype=np.uint8)\n",
    "mask2[pts[:, 0], pts[:, 1]] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill the right side this time, again using `skimage.measure.label` - the label of interest is 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask2 = flood_fill(mask2, (rmax, cmax), 1, connectivity=1)\n",
    "\n",
    "# The result\n",
    "plt.imshow(mask2, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final mask\n",
    "\n",
    "The last mask for the middle image is one of exclusion - it will be displayed everywhere `mask0` and `mask2` are not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask1 = ~(mask0.astype(np.bool) | mask2.astype(np.bool))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a convenience function to place masks in alpha channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_alpha(img, mask=None):\n",
    "    \"\"\"\n",
    "    Adds a masked alpha channel to an image.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    img : (M, N[, 3]) ndarray\n",
    "        Image data, should be rank-2 or rank-3 with RGB channels\n",
    "    mask : (M, N[, 3]) ndarray, optional\n",
    "        Mask to be applied. If None, the alpha channel is added\n",
    "        with full opacity assumed (1) at all locations.\n",
    "    \"\"\"\n",
    "    from skimage.color import gray2rgb\n",
    "    if mask is None:\n",
    "        mask = np.ones_like(img)\n",
    "        \n",
    "    if img.ndim == 2:\n",
    "        img = gray2rgb(img)\n",
    "    \n",
    "    return np.dstack((img, mask))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain final, alpha blended individual images and inspect them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pano0_final = add_alpha(pano0_warped, mask0)\n",
    "pano1_final = add_alpha(pano1_warped, mask1)\n",
    "pano2_final = add_alpha(pano2_warped, mask2)\n",
    "\n",
    "compare(pano0_final, pano1_final, pano2_final, figsize=(12, 12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have here is the world's most complicated and precisely-fitting jigsaw puzzle...\n",
    "\n",
    "Plot all three together and view the results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "# This is a perfect combination, but matplotlib's interpolation\n",
    "# makes it appear to have gaps. So we turn it off.\n",
    "ax.imshow(pano0_final, interpolation='none')\n",
    "ax.imshow(pano1_final, interpolation='none')\n",
    "ax.imshow(pano2_final, interpolation='none')\n",
    "\n",
    "fig.tight_layout()\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fantastic! Without the black borders, you'd never know this was composed of separate images!\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bonus round: now, in color!\n",
    "\n",
    "We converted to grayscale for ORB feature detection, back in the initial **preprocessing** steps. Since we stored our transforms and masks, adding color is straightforward!\n",
    "\n",
    "Transform the colored images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identical transforms as before, except\n",
    "#   * Operating on original color images\n",
    "#   * filling with cval=0 as we know the masks\n",
    "pano0_color = warp(pano_imgs[0], (model_robust01 + offset1).inverse, order=3,\n",
    "                   output_shape=output_shape, cval=0)\n",
    "\n",
    "pano1_color = warp(pano_imgs[1], offset1.inverse, order=3,\n",
    "                   output_shape=output_shape, cval=0)\n",
    "\n",
    "pano2_color = warp(pano_imgs[2], (model_robust12 + offset1).inverse, order=3,\n",
    "                   output_shape=output_shape, cval=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the custom alpha channel masks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pano0_final = add_alpha(pano0_color, mask0)\n",
    "pano1_final = add_alpha(pano1_color, mask1)\n",
    "pano2_final = add_alpha(pano2_color, mask2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View the result!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "# Turn off matplotlib's interpolation\n",
    "ax.imshow(pano0_final, interpolation='none')\n",
    "ax.imshow(pano1_final, interpolation='none')\n",
    "ax.imshow(pano2_final, interpolation='none')\n",
    "\n",
    "fig.tight_layout()\n",
    "ax.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the combined, color panorama locally as `'./pano-advanced-output.png'`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.color import gray2rgb\n",
    "\n",
    "# Start with empty image\n",
    "pano_combined = np.zeros_like(pano0_color)\n",
    "\n",
    "# Place the masked portion of each image into the array\n",
    "# masks are 2d, they need to be (M, N, 3) to match the color images\n",
    "pano_combined += pano0_color * gray2rgb(mask0)\n",
    "pano_combined += pano1_color * gray2rgb(mask1)\n",
    "pano_combined += pano2_color * gray2rgb(mask2)\n",
    "\n",
    "\n",
    "# Save the output - precision loss warning is expected\n",
    "# moving from floating point -> uint8\n",
    "io.imsave('./pano-advanced-output.png', pano_combined)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "<div style=\"height: 400px;\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"height: 400px;\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Once more, from the top\n",
    "\n",
    "I hear what you're saying. \"Those were too easy! The panoramas had too much overlap! Does this still work in the real world?\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Go back to the top. Under \"Load Data\" replace the string `'data/JDW_03*'` with `'data/JDW_9*'`, and re-run all of the cells in order.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "<div style=\"height: 400px;\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext load_style\n",
    "%load_style ../../themes/tutorial.css"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
